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Abstract: In this paper, we present a blind image restoration algorithm to reconstruct a high 
resolution (HR) color image from multiple, low resolution (LR), degraded and noisy images 
captured by thin (< 1mm) TOMBO imaging systems. The proposed algorithm is an ex- 
tension of our grayscale algorithm reported in [1] to the case of color images. In this color 
extension, each Point Spread Function (PSF) of each captured image is assumed to be differ- 
ent from one color component to another and from one imaging unit to the other. For the task 
of image restoration, we use all spectral information in each captured image to restore each 
output pixel in the reconstructed HR image, i.e., we use the most efficient global category 
of point operations. First, the composite RGB color components of each captured image are 
extracted. A blind estimation technique is then applied to estimate the spectra of each color 
component and its associated blurring PSF. The estimation process is formed in a way that 
minimizes significantly the interchannel cross-correlations and additive noise. The estimated 
PSFs together with advanced interpolation techniques are then combined to compensate for 
blur and reconstruct a HR color image of the original scene. Finally, a histogram normaliza- 
tion process adjusts the balance between image color components, brightness and contrast. 
Simulated and experimental results reveal that the proposed algorithm is capable of restor- 
ing HR color images from degraded, LR and noisy observations even at low Signal-to-Noise 
Energy ratios (SNERs). The proposed algorithm uses FFT and only two fundamental image 
restoration constraints, making it suitable for silicon integration with the TOMBO imager. 
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1. Introduction 

TOMBO (Thin Observation Module by Bound Optics) imaging systems are a new generation of 
paper-thin imaging products, which integrate micro-optics, photo-sensing elements and processing- 
circuitry all on one single silicon chip (Figure 1). These paper-thin imagers exhibit a thickness of less 
than a millimeter because they do not image the scene through a single imaging lens but through an array 
of microlenses [2]. The concept of a TOMBO imager was proposed and demonstrated by Tanida et al. 
[3-8]. The structure of a TOMBO imager is shown in Figure 1. It consists of an array of imaging units, 
each comprises a microlens that is associated with a subset of photo- sensitive pixel array. Individual 
imaging units are optically isolated by an opaque wall to prevent crosstalk (Figure 1). As a result, each 
individual imaging unit visualizes part of the scene. The output of the TOMBO imager is thus a mosaic 
of low resolution (LR) unit images. To reconstruct a high resolution (HR) image from the acquired set 
of LR images, Tanida et al. first proposed a Back-Projection (BP) method [6], which requires complete 
knowledge of the imaging system point spread function (PSF). The HR image of the scene is obtained 
by multiplying the captured LR images by the inverse (pseudo-inverse) of the known PSF. Tanida et al. 
proposed a second image reconstruction method (the "pixel rearrange method") [7], which computes 
the cross-correlation peaks between captured unit images to arrange and align unit image pixels in the 
HR image of the scene. A de-shading pre-processing step is introduced to compensate for the shading 
introduced by the separation walls (Figure 1). 



Figure 1. The architecture of a color TOMBO imaging system. 




We have previously proposed a novel spectral-based image restoration algorithm that require neither 
prior information about the imaging system nor the original scene [1]. Furthermore, the proposed al- 
gorithm alleviates the need for de-shading and rearrangement processing. In this paper, we extend this 
algorithm to color images. We examine the difference in characteristics between grayscale and color 
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images to develop a model for the color TOMBO imager. Previous work on color TOMBO imagers 
directly applied grayscale HR reconstruction algorithms to color images without considering the cross- 
correlation between color channels, and thus resulted in color artifacts [9-13]. In this paper, we exploit 
the global category of point operations for image restoration (see Figure 2) [14]. In this process, each 
pixel of the restored image is obtained by using information (pixels) from all captured LR images [15- 
20]. 



Figure 2. Point operations categories. 




Operations based on object structure/shape Geometrical Operations 




The proposed spectral-based color image restoration method averages out all LR captured images, 
making the color channels globally independent of each other. Compared to previously reported color 
restoration techniques [9], this proposed algorithm uses FFT and only two fundamental image restora- 
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tion constraints, which makes it suitable for silicon integration with a TOMBO imager. The paper is 
organized as follows: in section 2, we develop a model for a color TOMBO imager and derive the math- 
ematical formulation of a novel blind color image restoration method. Section 3 details the color image 
restoration algorithm. Results are discussed in section 4. Finally, a conclusion is given in section 5. 

2. Image Restoration Method for TOMBO Color Imaging Systems 

In this section, we extend the grayscale image restoration algorithm reported in [1] to color TOMBO 
imagers. In the restoration process, we consider the global point operations based on multiple images. 
By using this category of point operations, we exploit all available information in the mosaic of simul- 
taneously captured color images (see Figure 2). In addition, the global category is reported to have the 
ability to remove significant additive noise [15-20]. 

2. 1. System Model 

Consider a TOMBO color system with (/x x /j) captured color images as shown in Figure 1. Each 
captured color image represents a blurred, LR and noisy observation of an original HR scene. The 
mathematical model (Figure 3) for the system can be described by 
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(1) 

• gij(x, y, •&), ■d G {R, G, B] represents the blurred, LR and noisy color component for the i th J th 
captured unit image with resolution (M x N) pixels per color 

• hij(x,y,'d) is an (£ x £) PSF that represents the overall channel blur affecting g i! j(x,y,'d) unit 
image for the color component also called the intrachannel. We assume here that the blur is 
different for each color of each unit image 

• h ij j(x, y, GR), hij(x, y, BR), hij(x, y, BG) axe (£ x I) PSFs representing the overall mutual re- 
lation between red-green, red-blue and green-blue respectively. 

• " * *" represents the 2-D convolution operator w.r.t x, y 

• f(x, y, •&) is the t9 color component of the original scene with resolution (M x AT) > (M x N) 
per color component 

• Vi t j(x, y, "&) is the additive 2-D, zero mean white Gaussian noise per color component that affect 
the unit image g itj (x, y, ■&) 

• I D is the down-sampling operator representing the LR process 

Our main goal is to develop a restoration method that is able to reconstruct a HR, noiseless, color im- 
age of the original scene using only the (fj, x /j) LR, blurred and noisy TOMBO color images g>ij(x, y, 
The proposed method will have the following characteristics: (i) it does not require prior information 
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about the imaging system nor the original scene, and thus performs blind image restoration, (ii) it can re- 
move blur and additive noise from the HR color image, (iii) it exploits all available information contained 
in the captured LR images, and (iv) it requires minimal computational complexity. 



Figure 3. System model for the color TOMBO system. 

vi,i(x, y, i?) + cross-talks 




2.2. Formulation of the Restoration Method 

The general model of the color TOMBO imaging system is described by Equation (2). Let us consider 
all signal terms involved in the captured red color component of a unit image 

g ld (x,y,R) = [hij(x,y,R)**f(x,y,R) 

+ h itj (x,y,GR) **f(x,y,G) 

+ hi d (x,y,BR)**f(x,y,B) 

+ v id (x,y,R)] i D, i,j = l,..,fi 

By applying the two dimensional Z-transform to the model in Equation (2), we have 

Gij(zi,z 2 ,R) = £ XEHij(z£e^,z?e^,R)F(z*e^,z*e^,R) 

k=0 i=o v 7 v 7 

+ E H l Jz{>e^,zi>e^,GR)F(zfe^,z 2 D e^ J GR) 

k=0 1=0 v 7 v 7 /g\ 

+ E ZH id (zfe^,zfe^,BR)F(z?e^,z?e^,BR) 
k=o i=o v 7 v 7 

k=0 1=0 v 7 

where z\ = e j27r ^, ^2 = e^ 2n ^ 2 andGij(zi, z 2 , R), is the Z-transform of gi t j(x, y, R). 

The system model in Equation (3) is partitioned into the following terms: 

G id (z u z 27 R) = T^fz^z^R^+TP^z^z^R^+Tl^z^z^Rj 
+ Tfj (*? , zf , GR) + Tfj [zf , z? , BR) 
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where, 



T tj ( z f > z ? > = F *,i (V > ^ , ^ > ^ > + Hi , ^ , R\ represents the image of 
interest plus the noise term (defined in-frequency band useful terms), 

^(ztzhR) = D ± l EH iJ (zfe^,zpe^,R)F( K zre^,zh^,R),ar Q 1iio 
aliasing out of frequency band image terms, 



T?j {zf,z 2 D ,Rj = E E Vij ( zf e~~D— ? z£ e~ °~ ; , are the aliasing out-of-frequency band 



noise terms. 

l l \ fl-lD-l 



Tf tj (z*,z?,GR) = E E #y (zfe^^fe^iGR) F [zf e^, zfe^ 1 , d), are the 



fc=l !=1 

GR overall cross-talk terms. 

i i \ D-XD-l 



T^{zf \z° \BR) = E ±H i j(z?e^,z°e^,BR)F( K zTe^,z°e^,B), are 



fe=i i=i 

the BR overall cross-talk terms. 



In the above terms, the constant jp is not shown for simplification. 

By multiplying both sides of Equation (4) with the complex conjugate of the red blurring PSF, 

Hij [zf , zf , Rj , i.e., H*j (zf , zf , Rj and similarly by F* (zf , zf , Rj and applying the ensemble 
average operator, E {}, we have, 

e {f /, ,//;., (zp , , R) } = e {'/;;,//;., (zp , , i?) } + e {'/;':,//:, (zp , zp , r) } + E { /;.,//,;, (zf , zp , i? 

+ e {TfjHij hp , z* , gr) } + e {/;.,,///., (zp , zf , Si?) } 

E [G itj F* (zp , zp ,R)) = E {l^.F* (zp ,zp,R)}+E [t^F* (zp , zp , r) } + E {l^F* (zf , zf , r) } 

+ E [T^F* (zp , zp , Gr) } + E {l^ -F* (zp , zp , Br) } 



(5) 



which leads to: 



O ;..//• (zf , zf , R) = F (zf , zf , R) C Hi . Ht (zf , zf , r) + C VijH * (zf , zf , r) 

+ Qr* . Hlj (zf , zf , Rj + C Tli u h (zf , zf , Rj 

+ C Tfj Htj (zp,zp, GR) + C Ttj H h (zp, zp.BR) 

C Gi . F . (zp , zp , R) = Hl 3 (zp , zp , Rj C FF , (zp , zp , r) + C v ^f> (zf , zf , 



(6) 

C T * jF , (zp , z* , Rj + Ct^f* (zf , zp , i?) 
+ C T ^ F , (zf , zp , GR) + C T e. F , (zp , zp , BR) 

where cross-spectra Cxy*(z\, z 2 , R) = E {X(z±, z 2 , R)Y*(zx,z 2 , R)} over z\ y z 2 . 

Since the frequency bandwidth of the blurring PSF function H it j(zi,z 2 , R) is limited [15, 19, 21], the 
terms (zf , z 2 D > Rj an d (^i D > z 2 > ^ are not m tne same frequency band where JT*^ , 2 2 D , Rj 
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is. Therefore, the last two cross-spectral terms of Equation (6) can be discarded, leading to: 

C<;. //; , [z? ,z£,R) = F (zf , zf , i?) C n „■ (zf , zf , r) + C Vi>jHlj (zf ,z? ,R 

+ C T * j Ht . (zf , z * , Gi?) + C T *. h i i (zf , z* , BR} 

C GitjF . (zf , 4 > A) = flij (4 > . Cff* {4 , 4 > fl) + CV w f • (zf , z* , i?) 

+ Cjhi . p.. (zf , z * , Gi?) + Ct^ .f* (z? , 4 , BR} 



(7) 



i 



To minimize the effects due to interchannel cross -correlation terms (cross-talks) C T d H * ( zP , z 2 D ,R 

and Ct? h* ( z \ i z 2 > R ) m Equation (7), a typical approach would be to apply a decorrelation stretch 
technique to the captured LR images [14, 22]. However, such technique uses principal component anal- 
ysis, which is computationally expensive and thus not suited for silicon integration with a TOMBO 
imager. An alternative solution is to minimize the cross-correlation terms by using ergodicity, one of the 
wide-sense stationary properties of signals [23]. Let us consider the cross-correlations between the color 
channels, which depend on the reflection characteristics of the imaged objects. If the object exhibits 
strong reflection at the region of correlated wavelengths, then the cross-correlation between the channels 
could be strong. However, since the characteristics of one object is almost independent of those of other 
objects, the averaged interchannel cross-correlations over the entire image can be very weak. As a result, 
the RGB components of a color image can be regarded as being locally correlated but globally indepen- 
dent of each other [23]. In this paper, we exploit this property by using a global category of point opera- 
tions for image restoration. In this process, each pixel of the restored image is obtained by using informa- 
tion (pixels) from all captured LR images (Figure 2). This is carried out by computing the cross-spectra 
between the different color channels. Since the cross-spectra are nothing but the Fourier transform of 
the cross-correlations, our algorithm essentially averages out the interchannel cross-correlations over 
the entire image. As a result, the averaged interchannel cross-correlations become very weak [23]. Our 
spectral-based restoration method has the effect of further minimizing the interchannel cross-correlations 
because it averages out over all captured LR images. As a result, the cross-spectral terms of Equation 
7, C T a. Hh [zf , zf , GR^j the C T e. H *.* (zf , zf , BR) , C T a. F , (zf , zf , GR^j , C T ?. F * (zf , zf , BR 
can be neglected, leading to: 



C<;,ji-. (zP ,z 2 D ,R)=F (zP , zP , R) C Hi . Hl . (zP , zP , r) + C v . jHlj (zP ,zP,R 
GdjF* (zP , Z<? , R\ = Hij (zP , zP , Rj C ff * (zP , zP , Rj + Cy ijF * (zP ,zP,R 



(8) 



where Cv itj H^ (zP , zP , B^j an d Gy i jF * (z? , zP , R^j are the residual errors representing the in-band 
cross spectral terms between original images, in the PSFs and in the independent additive noise respec- 
tively. These errors can be assumed 2-D, independent and identically distributed (i.i.d.) signals under 
some regularity conditions [24]. By extending the previous analysis to the two other color components 
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and similarly 
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2.3. Restoration Process 



From Equations (8), (9), and (10), if some prior information or constraints about either the original 
or the blur is given, it is then possible to restore the original image. This can be done independently for 
each color component or jointly for all color components using Equations (9), and (10). For example, if 



the PSF associated with each captured g it j (x, y, R) unit image, H it j iz^ , z 2 D , R ) , is pre-determined or 
can be estimated, then it becomes possible to restore the original image using, 



F[zf,z 2 °,R 



z l i z 2 i 



R 



C HuH * [z°,z°,R 



(11) 



The above equation is valid under some constraints [24, 25]. One constraint is that the number of 
zeros in Hij(zi, z 2 , R) is much less than the size of the image. Because our method sums all PSFs 
spectra when implementing Equation (11), this will minimize the probability of having a division by 
zero. In general, if zeros exist but only in small numbers compared to the image size, a tolerance value 
can be added to the denominator to avoid division by zero. Other problems associated with image and 
PSF restoration methods are (i) the length £ and the PSF are unknown and most critically (ii) the impact of 
residual error terms, which can significantly deteriorate the restoration process [26, 27]. However, since 
a TOMBO imaging system can capture multiple observations of the same scene, it is possible to reduce 
the effects of the error terms significantly by using averaged cross-spectral techniques [25]. The use of 
global point operations will also minimize the additive noise [15-20]. To clarify this point, consider 
an imaging system with (// x //) captured images, the averaged spectral and cross-spectral techniques 
can be then applied to provide results similar to Equation (8) using spectral estimates instead of the true 
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estimates. In mathematical form, F (zf , z 2 D , R^j and similarly {^z° , z 2 D > ^ can be estimated using 
the averaged cross-spectra defined by 

A 1 A 4 / 1 1 \ M M a /ii 



i=lj=l v 7 i=l,j=l v 7 



where C*xy* (zi, -R) = ^(21,22, z 2 , -K) is an estimate of the cross-spectra between X(z 1? z 2 , 

and Y(z\, z%, R). For sufficient number of \i images, the last summation (error term of Equation (12)), is 
nothing but the mean value of an i.i.d. signal which has a zero mean [24, 25]. In this situation and after 
employing interpolation techniques, we have: 

A* A 1 „ 

E S L C G i , j Hl j {z 1 ,Z 2 ,R) 

F (21, ^2, i2) « : (13) 

EE^fl,', (zi,z 2 ,R) 

i=l j=l 

Under some constraints [17, 19, 20, 22, 25-28], the above estimates can be used to restore the original 
image. 

3. Color Image Restoration Algorithm 

In this section, we employ the work formulated above to develop a color image restoration algo- 
rithm. Using Equation (13), it is possible to restore the original image using iterative techniques. During 
the restoration process, the algorithm will only impose two fundamental image restoration constraints 
(positivity and support region): 

• For restoring the image 

\f(x,y,R)\, <x,y>e [LM, LN] 

./(.<•■ </./?) { F M , \f(x,y,R)\>F 00 ,<x,y>e [LM,LN] (14) 

0 otherwise 

where, Fm is the mean value of a mesh of pixels in the region < x, y > 6 [LM, LN] surrounding 
symmetrically the pixel index (x,y), Foo is a threshold of a large value pre-defined for the case 
when the summation of the PSF spectra in Equation 13 is close to zero. 

For estimating the PSFs 

k j (x, y ,R) = l^ X ' V > R)l > <X ' V>e M (15) 

0 otherwise 

where, L > D is the interpolation or up-sampling factor needed to restore the HR image, (LM x LN) 
is the size of the restored image which can be greater or equal to the size of the original image (M. x J\f), 
and (£ x t) is the size of the estimated PSF. 

Practical considerations for the iterative algorithm implementation can be summarized as follows: 
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• Pixel amplitudes that reach values greater than 255 are scaled using the following histogram nor- 
malization, 

(/O, V, R) - fmin,R) 

fscaied(x, y,R) = ax — — — h b (16) 

\Jmax,R Jmin,R) 

where a and b are usually 255 and 0 respectively (but other values can be also used to adjust 
contrast and brightness), f ma x,R an d f m in,R are the maximum and minimum pixel values of the 
color component R. For improved performance, f m i nj R and frn a x,R are usually chosen as the 5% 
and 95% levels in the histogram distribution respectively (confidence interval). 

• The mean value of the input image(s) and the output image is to be maintained (note that there are 
twice as many green pixels as red/blue pixel for the Bayer filter) 

• To resolve the problem of having zeros or nulls in the spectra, the following equation for the 
interpolated f(x, y, R) is used: 

T l T l c G i , j HT tj (zi,Z2,R)-a 1 
F z 2 , R) » ' (17) 

£ E Ch^h*. (zi, z 2 , R) + a 2 
i=ij=i 

followed by: 

( Zl ,z 2 , R) = pp™ (z h z 2 , R) + (l- /3)F fa, z 2 , R) (18) 

where cti,a 2 are two positive numbers representing the extent of additive noise in the residual 
terms, (3 is a recursive stability factor controls the amount of information needed from posterior 
image estimates (0 < /3 < 1 ), n + 1 is the current iteration number. Typical values of (3 range 
between 0.5 to 0.9. Like in adaptive systems, the value of (3 can be adjusted to avoid such that 
impulsive-like outputs. 

• For initialization, one of the images is used as an initial estimate of the HR image. The up- 
sampling and interpolation process is done by zero-padding in the spatial domain between the 
image samples. Afterwards the FFT is applied. In the Fourier domain, a single spectrum is then 
taken out of the repetitive spectra using a low pass filter with cut-off frequency (j) and zeroing 
the rest of the spectrum. Finally, inverse fast Fourier transform (IFFT) is applied to inverse back 
to the image domain. It is essential that the zero-padding be done such that the zero frequency 
components remain the same. In addition, zero-padding should be applied to both positive and 
negative frequencies. Unlike existing techniques that use lower order functions for interpolation 
(cubic interpolation used in [2]), our method uses the more efficient sine function. 

• We use the 2-D fast fourier transform (FFT) to estimate spectra and cross spectra needed for the 
algorithm 

The proposed HR color image restoration algorithm is detailed in Table 1 . 
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3.1. Convergence Analysis 

From Equation (18) and after considering the situation where n = 0, 1, 2, . . . , n a (note that F also 
changes at each iteration), after some mathematical manipulation it can be shown that 

n a 

( Zl , z 2 , R) = p^+VpW (zi, Z2, R) + J2 W ~ P)F {n °~ j) (*i> z ^ R) (19) 

3=0 

For a value of 0 < (3 < 1 and when the number of iteration n Q — ► oo or large enough, the above 
equation can be approximated to: 

( Zl , z 2 , R) = J2 Pi 1 ~ P)F {n °~ j) R) (20) 

3=0 

Based on spectral estimation principals [24, 25], and by using the image estimator provided in Equa- 
tion (17), it can be proven that: 

e{f(z 1 ,z 2 ,R)}=F(z 1 ,z 2 ,R) (21) 

Thus 

E |#(n 0+ i) (^ 2)jR )} = J^Fil-P) e{f^ (22) 

i=o 

E{#("o+ 1 )( Zl , Z 2,i2)} = (l-^)F( Zl ,z 2 ,i2)E^ 
k J i=o 

= F (zi, z 2 , R) (1 — J (23) 

= F(z u z 2 ,R) 

From the above analysis, it can be clearly seen that for any value of 0 < (3 < 1, the algorithm will 
converge to an unbiased estimator of the original image. 

4. Results and Discussion 

To evaluate the performance of the proposed color image restoration algorithm, we consider the fol- 
lowing tasks: 

(1) Restore a HR image from multiple blurred, LR and noisy "simulated" TOMBO color images. 

(2) Restore a HR image from multiple blurred, LR and noisy images and compare the results with 
the previous method [9, 10]. 

(3) Restore a HR image from multiple blurred, noisy "real" TOMBO color images. 
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Table 1. Proposed color image restoration algorithm. 

Step 1: Set the values of L, £, a±, a 2 , (3, F M , 

Step 2: Select the color to be restored d G {R, G, B} 

Step 3: Iteration n = 1, initialize f(x, y, hence F(f 1 , f 2 , i?) = FFT y, #)j 
Step 4: For i, j = 1,2,...,// estimate 

C GtjF .(f 1 ,f2,#) = G i j(f 1 ,f 2 ,#)*F*(f 1 ,f 2 ,0) (24) 



HiAfi>M= NT^B^J ^■(x,y^)=IFFr{ff iJ -(/ 1 ,/ 2 ^)} (25) 

F(f 1 ,f 2 ,$)F*(f 1 ,f 2 ,$) + a 2 J 

Step 5: Impose PSF constraints to get the accurate estimates 

hij(x,y,#) = < Q otherwise j = I ^' / 

(26) 

Step 6: Estimate the biased image spectra 



EE^(/i,/ 2 ,*)-ai 
^ / 2 , 0) = =► / (x, y, 0) = IFFT {F (/i, / 2 , 0) } (27) 



i=ij=i 

Step 7: Impose the image constraints 

f |/>,y,tf)|, <x,y>G [LM, LN] 

f(x,y,&)=l F M , \f(x,y,0)\>F oo ,<x,y>e [LM,LN] (28) 
[ 0 otherwise 

then estimate the original image by updating the estimates using 

/ (x, y, 0) = 0f (x, y, ■&) + (1 - 0)f (x, y, ■&) (29) 

Step 8: Scale the estimated images pixels or use histogram normalization to find a and b, then adjust the 
image using 

(f(x,y,^)-fmin^) 

f(x, y, #) = a x ^— : ^ + 6 (30) 

( fmax,tf fmin;d ) 

Step 9: Repeat from Step 4 until convergence, then repeat for another color d 
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4. 1. Examples of Simulated Images 

In this section, we test the performance of our proposed algorithm in restoring a HR image from 
simulated TOMBO images of "Lena" [29]. The algorithm performance is compared with the pixel 
rearrange method developed in [2] for TOMBO imagers. Simulation parameters for the generated LR, 
blurred and noisy images are given in Table 2 for noiseless and noisy cases. The simulated images 
are generated following the procedure described in [1]. In all simulations, the SNER is defined as 



Table 2. Input parameters for simulated TOMBO images. 
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Figure 5 
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Figure 4 shows that our restoration algorithm can restore a HR color image from the simulated LR, 
blurred TOMBO color images in the absence of additive noise. Our method performs better than the 
pixel rearrange method because the pixel rearrange method can not align captured images. In Figure 5, 
we tested the algorithm with the simulation parameters given in Table 2, but with additive noise. Our 
algorithm appears almost insensitive to additive noise, whereas a significant amount of noise can still be 
seen in the image restored by of the pixel rearrange method. From the two simulation results, we can see 
that our proposed algorithm converges rapidly. 

4.2. Comparison with Existing Image Restoration Methods 

In this section, we compare our algorithm with the advanced restoration methods developed by Sina 
in [9] and Sroubek in [10]. We generated 12 blurred, LR, noiseless and noisy images following the 
procedure explained in [9]. Results for the two cases are shown in Figure 6, and the PSNR [9] for each 
restoration is summarized in Table 3. 



Table 3. PSNR values (in dB) for the restoration methods. 





Sina [9] 


Sroubek [10] 


This Work 


Noiseless 


21.986 


18.250 


16.348 


Noisy 
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Figure 4. Simulation results, 6x6 images, D = 4, L = 4, £ = 7, no additive noise. 
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Figure 5. Simulation results, 6x6 images, D = 4, L = 4, i = 7, SNER = 4.968 dB. 
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Figure 6. Simulation results for 12 blurred LR images of the lighthouse, in noiseless, and 
noisy, SNER=10.4616 dB, D = 4, L = 4, i = 5. 
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At high SNRs, our algorithm does not perform as well as Sina's, which enforces more constraints 
including (i) a penalty term to enforce similarities between the raw data and the HR estimate (data fi- 
delity penalty term), (ii) a penalty term to encourage sharp edges in the luminance component of the HR 
image (spatial luminance penalty term), (iii) a penalty term to encourage smoothness in the chrominance 
component of the HR image (spatial chrominance penalty term), and (iv) A penalty term to encourage ho- 
mogeneity of the edge location and orientation in different color bands (inter-color dependencies penalty 
term). In Sroubek's method, the following constraints are enforced: (i) a fidelity term, (ii) a smooth- 
ing term using variational integrals, (iii) a consistency term that binds the different volatile PSFs, (iv) a 
smoothing term to overcome the higher nullity of integer-factors, and (v) an anisotropic term for edge 
preservation. Although our algorithm only imposes two fundamental constraints, its performance is vi- 
sually satisfactory, as seen in Figure 6. Our method achieved a PSNR of 16.348 dB in contrast to a PSNR 
of 21.986 dB and 18.250 dB for Sina's and Sroubek's methods, respectively. A higher PSNR, however, 
does not always correlate to subjective quality [30, 31]. For instance, the Shift- And- Add restored image 
(in Sina et al. [9] page 148) has a PSNR of 17.17 dB, but is clearly of poor quality compared to the one 
we restored. At low SNRs (10.46 dB), our algorithm visually outperforms both aforementioned methods 
(Figure 6). However, their PSNRs (14.13 dB for Sina's and 13.78 dB for Sroubek's) are higher than the 
10.89 dB achieved by our algorithm. This can be explained by the fact that both Sina's and Sroubek's 
methods minimize regularized energy functions. 

Furthermore, we have observed that Sina's method experience some limitations when applied to im- 
ages blurred using semi-gaussian PSFs. This could be due to the gaussian kernel used in Sina's approach. 
It is also important to point out that Sina's method uses multiple frames of an image with different dis- 
placements (i.e. diversity is guaranteed, see section VI in [9]), while our method uses multiple simultane- 
ous observations of the same scene captured by the TOMBO imager. The diversity between all captured 
images cannot be guaranteed for the case of our TOMBO imager. 

4.3. Examples of Real Images 

In this section, we investigate the performance of our proposed algorithm with real captured TOMBO 
color images. In this example, the captured images are that of a "teddy bear" picture, with a plan object 
located at 200 mm from the sensor module. Each unit imager has 60 x 60 pixels and each pixel is 
6.25 [Lm x 6.25 [mi. The microlens array has the following characteristics: 1.3 mm focal length, 0.5 
mm diameter of aperture, and a 0.5 mm pitch for the microlens array. The TOMBO imager integrates 
a color filter array with a Bayer pattern. Demosaicing is achieved on-chip. To test the performance of 
the restoration algorithms at lower SNRs, a zero mean white Gaussian noise is added manually to the 
captured LR images. System parameters for this real example are given in Table 4. 



Table 4. Input parameters for real captured TOMBO images, D = 4. 
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Figure 7. Experimental results, 4x4 unit images, D = 4, L = 4, t = 5, SNER = 15.774 dB. 
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Figure 7 demonstrates that our proposed algorithm is able to restore a HR color image of the "teddy 
bear" picture using real captured TOMBO LR observations. In addition, our algorithm can significantly 
minimize the additive noise. It outperforms the conventional pixel rearrange method [33, 34]. Further- 
more, our algorithm can restore the HR color image within only 15 iterations. 

5. Conclusions 

A blind color image restoration method is proposed for the reconstruction of HR color images using 
multiple LR, degraded and noisy color images captured by TOMBO imagers. The proposed spectral- 
based method only imposes two fundamental image restoration constraints (positivity and support re- 
gion) to (i) correct the blur that affects the captured LR images, (ii) minimize the interchannel cross- 
correlations between RGB color components, (iii) significantly reduce the impact of additive noise, and 
(iv) reconstruct a HR color image. The computation complexity of the algorithm is low compared with 
existing techniques because it uses FFT for spectral estimation. 

The proposed restoration algorithm has a rapid convergence rate of 10 to 20 iterations. Results show 
that the proposed algorithm is capable of restoring a HR image from the degraded LR color TOMBO 
images even when the SNER is as low as 5 dB. The proposed algorithm uses FFT and only two funda- 
mental image restoration constraints, which makes it suitable for silicon integration with the TOMBO 
imager. 
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